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We extend our recently developed interatomic potentials for UO2 to the mixed oxide fuel sys- 
tem (U,Pu,Np)02. We do so by fitting against an extensive database of ah initio results as well 
as to experimental measurements. The applicability of these interactions to a variety of mixed 
environments beyond the fitting domain is also assessed. The employed formalism makes these 
potentials applicable across all interatomic distances without the need for any ambiguous splining 
to the well-established short-range Ziegler-Biersack-Littmark universal pair potential. We therefore 
expect these to be reliable potentials for carrying out damage simulations (and Molecular Dynam- 
ics simulations in general) in nuclear fuels of varying compositions for all relevant atomic collision 
energies. 



The interest in using Mixed Oxide (MOX) fuels com- 
prising (U,Pu,MA)02 (where MA = Np, Am and Cm) in 
fast breeder and transmutation reactors is ever increas- 
ing. Since this complex fuel experiences a high burn- 
up ratio with large quantities of fission products and 
materials defects, it becomes crucial to understand the 
evolution and statistics of atomic displacement cascades 
due to high energy radiation that the material faces^. 
Classical Molecular Dynamics (MD) with its ability to 
simulate fairly long length scales, though still retain- 
ing the fine atomic structure of the material, is ide- 
ally suited for such studies. However, the complexity 
of the interatomic interactions for radiation damage sim- 
ulations cannot be fully represented by simple classical 
forms due to the disparate scales of energies involved. 
Interactions corresponding to equilibrium conditions are 
traditionally found by fitting to a variety of thermo- 
dynamic data; while for description of the short-range 
behavior, the Ziegler-Biersack-Littmark (ZBL) universal 
pair potential^ developed in the 1980s is well-accepted. 
These two "pieces" then need to be smoothly connected 
via somewhat arbitrarily applied splines. We recently 
proposed a methodology for developing interatomic po- 
tentials that is valid for all interatomic separations, with- 
out the need for any ambiguous splines^. In this article, 
we apply this formalism to a more general case of MOX 
nuclear fuels of varying composition. In addition to cap- 
turing high temperature thermodynamic properties, as 
done by available potentials^"'', we also incorporate cor- 
rect treatment of point defects. Created due to irradia- 
tion, these are critical for the understanding of a variety 
of phenomena such as fuel swelling, fission gas release and 
burn-up structure formatiortii^. A key test of any devel- 
oped energy surface lies in its ability to adequately rep- 
resent systems/configurations that were not included in 
the fitting procedure. - We here fit the potential param- 
eters to ab initio and experimental data for the oxides 
Pu02 and Np02, and then check for their transferabil- 
ity by comparing against ab initio data for (Ua;Pui_2,)02 
and (Ua;Npi_j,)02 configurations that were not included 



in the fit. 

In the present study we employ the generalized po- 
tential formalism- that behaves correctly in both short- 
range and long-range limits. The only component in 
this potential that remains to be determined is a cor- 
rection term for intermediate distances associated with 
chemical bonding. We find this correction term by fit- 
ting to an extensive database of generalized gradient ap- 
proximation GGA-I- U ab initio calculationsifi on Pu02 
and Np02. The potential's applicability in a mixed 
environment pertinent to MOX fuels is further veri- 
fied by testing against GGA+[/ data for (U2;Pui_a;)02 
and (Uj;Npi_a;)02. GGA+f/ is known to provide elec- 
tronic and magnetic behaviors of the actinide oxides^i 
that are consistent with experiments. In this approx- 
imation, the spin-polarized GGA potential is supple- 
mented by a Hubbard-type term to account for the 
localized and strongly correlated 5/ electrons. Our 
database comprises results obtained from GGA+ U cal- 
culations with the projector augmented-wave method 
and coUinear antiferromagnetic moments as implemented 
in the VASP packagei^,. Dudarev's rotationally invariant 
approacbi^ii^ to GGA-|- U is employed wherein the pa- 
rameter U-J is set to 3.99, 3.25 and 3.40 for U, Pu and Np 
respectivelyi^"— . These are the generally accepted values 
for reproducing the correct band structures of the corre- 
sponding oxides. Energy cutoff for the plane waves was 
kept at 400 eV. Since GGA-I- U overestimates the lattice 
parameter, a common scaling factor (same as that used^ 
for UO2) was employed to get experimentally correct lat- 
tice parameters. The ab initio database so obtained for 
fitting comprises: 

1. Isochoric relaxed runs on a 12 atom unit cell, which 
was isometrically contracted and expanded by var- 
ious amounts (i.e., equation of state calculations 
wherein each data point was calculated under the 
constraint of constant cell volume) and for which an 
8x8x8 fc-point grid was taken after ascertaining k- 
point convergence. Ionic relaxations were carried 
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out until residual forces were less than 0.01 eV/A. 

2. Static (i.e., no ionic relaxation) runs on a 96 atom 
2x2x2 supercell in which one atom at a time (O 
or Pu or Np) was perturbed from its equilibrium 
position by varying distances (on the order of 1 A 
or less from the equilibrium positions) in different 
directions. Sampling of the gamma point only was 
found to be satisfactorily accurate for this. 

3. A 96 atom 2x2x2 supercell for the formation en- 
ergies of stoichiometric defects, namely, Oxygen 
Frenkel pair. Neptunium Frenkel pair and Pluto- 
nium Frenkel pair. Several vacancy-interstitial dis- 
tances were considered to ascertain the separation 
between these corresponding to the minimum de- 
fect formation energy (excluding the case of near- 
est neighbor distances, which was found to lead to 
vacancy- interstitial recombination). Correct pre- 
diction of these energies has been given great im- 
portance in generating interatomic potentials for 
cascade simulations in UO^ii'^'iSri^. 

A total of approximately 50 ab initio configurations were 
thus used in the fitting. Note that in the above calcula- 
tions, any interactions between atoms and their periodic 
images do not systematically bias the fit of the potentials 
because the same supercell geometry is used in both the 
ab initio and the empirical potential energy calculations. 

The ab initio database employed for validation and 
for testing transferability includes equation of state runs 
similar to those in the fitting database, for oxides of 
UaiPu, U30PU2, UaiNp and U30NP2, each with 64 Oxy- 
gens. These data points were not included in the fit it- 
self and were used only after the fitting was complete for 
checking the robustness of the potentials with respect to 
use in mixed environments. 
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FIG. 1: Test of approximation- validity of fp„Pu = {90/88)fuu 
and fjvpjvp — {89/88)fuu by looking at the applicability of 
similar relations for cations of members of the previous row 
of the periodic table with similar shell structure viz. Pm and 
Nd. Dashed line denotes the result from this approximation 
while solid line is the actual charge density^ for Nd^** . 

In addition to the ab initio data, we also included ex- 
perimental thermal expansion behaviorS" of Pu02 and 
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FIG. 2: Quality of fit from our fitted potential (asterisks) 
for various ab initio energies (circles) for Pu02 :(a) equation 
of state (b) oxygen atom perturbation (c) plutonium atom 
perturbation. For each of oxygen and plutonium, the first 
four perturbations are along (100) direction while the second 
four are along (110) direction. The perturbations are on the 
order of 1 A or lower from the equilibrium positions. 
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FIG. 3: Quality of fit from our fitted potential (asterisks) 
for various ab initio energies (circles) for Np02 :(a) equation 
of state (b) oxygen atom perturbation (c) neptunium atom 
perturbation. For each of oxygen and neptunium, the first 
four perturbations are along (100) direction while the second 
four are along (110) direction. The perturbations are on the 
order of 1 A or lower from the equilibrium positions. 



Np02 in the fit. We found that including experimental 
thermal expansion data (which is readily available) is a 
very effective way to ensure reasonable thermal expan- 
sion behavior in this system. To make the calculation 
of high temperature lattice parameters computationally 
tractable during the fitting procedure, we employed the 
quasiharmonic approximation (QHA)^^, in which atoms 
are treated as pure harmonic oscillators whose frequen- 
cies depend on the cell volume. The so-called zero static 
internal stress approximation (ZSISA)^"^ to QHA, as im- 
plemented in the package GULP, was used^i. QHA in- 
volves a full relaxation with respect to external (cell pa- 
rameters) and internal (atom positions within the cell) 
coordinates. ZSISA ignores the dependence on internal 
coordinates of the vibrational part of the free energy. 
We found that for the materials studied and potential 
forms used in this communication, the lattice parame- 
ter through NPT (constant Number, Pressure, Tempera- 
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ture) MD was slightly lower than that through ZSISA. As 
such, an empirical adjustment to the ZSISA lattice pa- 
rameter had to be included in the fitting. Thus, several 
independent fits were done using ZSISA lattice parameter 
values equal to the experimental lattice parameter multi- 
plied by 1], with 77 varying between 1 and 1.01. NPT MD 
was carried out with these potentials (details of MD pro- 
vided later) to find the 77 that led to MD values matching 
the experimental data the best. We found that rj equals 
1.0006 and 1.0008 for PuOz and NpOa respectively, for 
a best match in the least squares sense between experi- 
mental and NPT MD lattice parameters. 
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FIG. 4: (a) fitted O-Pu interaction (b) fitted O-Np interaction 



The potential forms thus used for fitting to the ab ini- 
tio and experimental data are similar to that proposed 
previouslj*^, and are summarized below for Pu-Pu and 
Pu-0 interactions (with similar forms for other interac- 
tions): 



VoPu 



VpuPu{r) = ZBL,oMr)+^-^^ + ^[--—fpuPu{r)] V 0<r (1) 

r ZBL,oMr) + ^[f - ^fooir)] ~ ^[f - ^fPuPu{r)] < r < r, 

/\ _ ( ^)(4)e ^ I ^th Qj,(igj. polynomial ri < r < r2 

~ Aircor | Aexp{-r/ p) - B /r^ + {r - r3)'^{Cr^ + Dr^) r2 < r < 
[ Aexp(— r/p) — B/r^ < r 



TABLE I: Defect energy comparisons 





ab initio 


Potential 


Potential 




(Current work) 


(Current work) 


(Previous works^i^) 


O Frenkel pair formation energy in Pu02 (eV) 


3.9 


4.9 


7.0 


O Frenkel pair formation energy in Np02 (eV) 


4.5 


5.8 


10.0 


Pu Frenkel pair formation energy (eV) 


11.9 


24 


17 


Np Frenkel pair formation energy (eV) 


12.2 


26.7 


17.5 



The UO2 family of interactions is kept the same as in 
Rcf. 3. Here ZBLz, +gi_22-i-g2 (^) denotes the ZBL form 
of interaction between two neutral atoms having atomic 
numbers Zi + qi and Z2 + 52 , but using the screening 
length for Zi and Z2, as explained in Ref. i3|. The func- 



tions / in the above are related to the charge densities 
of the respective atoms. Detailed coefficients of ioo and 
f[/[/ can be found in Ref. U while fp„p„ and iNpNp can 
be calculated from the relations ipuPu = {Q0/88)iuu and 
iNpNp = (89/88)f£/(7. This was needed since Np+^ and 
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FIG. 5: Equation of state for (a) U31PUO64 and (b) 
U3oPu2 64. Circles denote ab initio data while asterisks are 
the values predicted (not fitted) with current potential. 



-933 p 
-933.5 

-934 
-934.S3 

-935 
-935.5 



(a) 



O ab initio 

¥ from current potential 



O 
* 



* 

o 



o 

* 



-0.5 0.5 

% change in cell parameter 



-932.6 p 
-932.9= 

-933 
-933.?) 
-933.4 
-933.6 
-933.8 

-934 
-934.2 



(b) 



O ab initio ( 
^ from current potential 



O 
* 



o 
* 



-0.5 0.5 

% change in cell parameter 



FIG. 6: Equation of state for (a) U3iNp064 and (b) 
U30NP2O64. Circles denote ab initio data while asterisks are 
the values predicted (not fitted) with current potential. 



Pu"*"^ charge densities p{r) are not available in Ref. 0. 
Wc tested this approximation using cations of elements 
in the previous row of the periodic table where actual 
ZBL charge densities are available, viz. Nd, Pm and Sm. 
As can be seen from figure [H the approximation satisfac- 
torily captures the electronic shell structure of 47rr^p(r), 
which is the quantity of interest to us. Note that we have 
removed altogether any splines for cation-cation interac- 
tions. The downhill simplex method of Nelder-Mead was 
then used to carry out the potential fitting^^. The fit- 
ting involved minimizing an objective function equal to 
the sum of the squares of the differences between the 
ab zmtio/experimental data (weighted since they denote 
different quantities) and that predicted by the poten- 
tial for all the classes of data as detailed above. GULP 
was used for energy calculations and for atomic-positions 
optimization^i. 

Figure [2] shows the quality of fit for the Pu02 equa- 
tion of states and single atom perturbation data, while 
figure [3] shows the same for Np02. Table [T] shows the de- 
fect formation energies as obtained by us in the GGA-I- U 
calculations, along with the corresponding values from 
the current potential and from the previous potentials 
published for these systems. We excluded the cation de- 
fect formation energies entirely from the fitting objective 
function. This can be justified by considering that (i) 
these energies as per ab initio are already very high - up- 
wards of 12 eV; (ii) from the case^ of UO2, it is expected 
that ab initio actually underestimates these energies, and 
thus they are even less likely to form; and (iii) these (Pu 
and Np) are the minority cations. It has been argued^^ 
though that Uranium Frenkel pairs and Schottky trios 
might play an important role in the diffusion of noble 
gas impurities formed after fission - as such, our library 
of potentials does provide a much better match for the 
Uranium Frenkel pair and Schottky trio formation energy 
since it is based on the potentials in Ref. 

The potentials so obtained are plotted in Figure SI 
while the fitted coefficients are detailed in Table HIl Note 
that since there was no spline in any cation-cation inter- 
action (see Equation H])), they do not find a mention in 
the above list. The aforementioned S*'' order polynomial 
is uniquely determined by the provided cutoffs and po- 
tentials. The detailed potentials are available as a GULP 
library file. 



TABLE II: Coefficients of fitted potentials 





PUO2 


NpOa 


A (eV) 


597.304 


597.605 


P(A) 


0.475712 


0.484948 


B (eVA") 


0.31187 


0.31187 


C (eV/A") 


0.0003375 


-0.0735556 


D (cV/A*) 


0.029085 


0.048972 


ri (A) 


1.42 


1.17 


r2 (A) 


1.7 


1.7 


ra (A) 


2.85 


2.94 
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FIG. 7: Lattice parameter at various temperatures for (a) 
Pu02 and (b) Np02. Straight lines are the experimental 
values^ valid between 400 and 1000 K, while circles denote 
values obtained from MD simulations using current poten- 
tials. Plus signs represent (l/rj) times the experimental val- 
ues actually used in fitting to account for the observation 
that ZSISA slightly overestimates the MD lattice parameters. 
Details of calculation of this adjustment factor rj (equaling 
1.0006 and 1.0008 for Pu02 and Np02 respectively) can be 
found in the text. 



and experiments is excellent. The quality of the enthalpy 
values compared between experiments^! and those pre- 
dicted from NPT MD with current potential is also very 
good (see figure 

To summarize, we have developed interatomic potcn- 
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FIG. 8: Enthalpy at various temperatures (relative to room 
temperature enthalpy) for (a) Pu02 and (b) Np02. The 
circles denote values from NPT MD (predicted and not fit- 
ted values) while the asterisks are the known experimental 
values^!. 



The performance of the potential against the validation 
data, i.e., equation of states for oxides of UsiPu, U30PU2, 
UaiNp and U30NP2 can be seen from figures [S] and HI The 
match is satisfactory and interestingly it improves with 
more Pu or Np content in respective cases. 

The generated potentials were verified through NPT 
MD simulations on 3x3x3 unit cells (324 ions). The 
system was equilibrated for 10 ps while production runs 
were carried out for 100 ps with time steps between 0.001 
and 0.0005 ps (depending on temperature). Apart from 
the lattice parameter, we also considered the enthalpy as 
a function of the temperature. 

Figure [7] compares the lattice parameter as obtained 
from the MD simulations with experimental values for 
Pu02 and Np02^^. Figure [7] also shows the correspond- 
ing ZSISA values as obtained from the potentials. The 
over-estimation adjustment factor 77 used on the ZSISA 
values can be seen here. After this adjustment to ZSISA, 
the match for the lattice parameters between NPT MD 



tials for the Mixed Oxide fuel system (U,Pu,Np)02 by 
fitting to an extensive ab initio database and to available 
experimental observations using a formalism that has 
been shown to be capable of dealing in a self-contained 
manner with conditions ranging from thermodynamic 
equilibrium to very high energy collisions relevant for fis- 
sion events. The potentials capture known experimental 
measurements on these oxides as well as a rich database 
of ab initio GGA-I- U results. The applicability of these 
potentials in scenarios not included in the fitting is also 
explicitly demonstrated. 

This research was supported by the US National Sci- 
ence Foundation through TeraGrid resources provided by 
NCSA under grant DMR050013N, through the U.S. De- 
partment of Energy, National Energy Research Initiative 
for Consortia (NERI-C) grant DE-FG07-07ID 14893, and 
through the Materials Design Institute, Los Alamos Na- 
tional Laboratory contract 75782-001-09. 



^ R. Devanathan, L. Van Brutzel, A. Chartier, C. Gueneau, and P. Van Ufi^elen, Energy Environ. Sci. 3, 1406 (2010). 

A. Mattson, V. Tikare, T. Bartel, T. Besmann, M. Stan, 



6 



^ J. F. Zicglcr, J. P. Biersack, and U. Littmark (Pergamon, 

New York, 1985). 
^ P. Tiwary, A. van dc Walle, and N. Gr0nbech-Jensen, Phys. 

Rev. B 80, 174302 (2009). 
■* K. Kurosaki, K. Yamada, M. Uno, S. Yamanaka, K. Yar 

mamoto, and T. Namekawa, Journal of Nuclear Materials 

294, 160 (2001). 
^ T. Arima, S. Yaniasaki, Y. Inagaki, and K. Idcmitsu, Jour- 
nal of Alloys and Compounds 400, 43 (2005). 
® K. Kurosaki, M. Imanmra, 1. Sato, T. Namekawa, M. Uno, 

and S. Yamanaka, Journal of Alloys and Compounds 387, 

9 (2005). 

N. D. Morelon, D. Ghalcb, J.-M. Delaye, and 

L. Van Brutzcl, Phil. Mag. 83, 1533 (2003). 
® L. Van Brutzel, A. Chartier, and J.-P. Crocombette, Phys. 

Rev. B 78, 024111 (2008). 
® A. van de Walle and G. Ceder, Journal of Phase Equilibria 

23, 348 (2002). 

V. I. Anisimov, J. Zaanen, and O. K. Andersen, Phys. Rev. 
B 44, 943 (1991). 
" K. T. Moore and G. van der Laan, Rev. Mod. Phys. 81, 
235 (2009). 

G. Kresse and J. Furthmiiller, Phys. Rev. B 54, 11169 
(1996). 

S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. 
Humphreys, and A. P. Sutton, Phys. Rev. B 57, 1505 
(1998). 

" S. L. Dudarev, G. A. Botton, S. Y. Savrasov, Z. Szotek, 



W. M. Temmerman, and A. P. Sutton, Phys. Stat. Sol. A 
166, 429 (1999). 

H. Y. Geng, Y. Chen, Y. Kaneta, M. Iwasawa, T. Ohnuma, 
and M. Kinoshita, Phys. Rev. B 77, 104120 (2008). 
H. Y. Geng, Y. Chen, Y. Kaneta, and M. Kinoshita, Phys. 
Rev. B 75, 054111 (2007). 

D. A. Andersson, J. Lezama, B. P. Uberuaga, C. Deo, and 

S. D. Conradson, Phys. Rev. B 79, 024110 (2009). 

R. Devanathan, J. Yu, and W. J. Weber, The Journal of 

Chemical Physics 130, 174502 (2009). 

K. Covers, S. Lemehov, M. Hou, and M. Verwerft, Journal 

of Nuclear Materials 366, 161 (2007). 

K. Covers, S. Lemehov, M. Hou, and M. Verwerft, Journal 

of Nuclear Materials 376, 66 (2008). 

T. Yamashita, N. Nitani, T. Tsuji, and H. Inagaki, Journal 
of Nuclear Materials 245, 72 (1997). 

22 A. van de Walle and G. Ceder, Rev. Mod. Phys. 74, 11 
(2002). 

23 N. L. Allan, T. H. K. Barron, and J. A. O. Bruno, The 

Journal of Chemical Physics 105, 8300 (1996). 
2" J. D. Gale, J. Chcm, Soc, Trans. 93, 629 (1997). 
"^'^ J. A. Nelder and R. Mead, The Comp. Jour. 7, 308 (1965). 
^® A. Chartier, L. Van Brutzel, and M. Freyss, Phys. Rev. B 

81, 174111 (2010). 

H. Serizawa, Y. Aral, and K. Nakajima, The Journal of 
Chemical Thermodynamics 33, 615 (2001). 



